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Abstract 

The osmotic virial coefficient B2 of globular protein solutions is calculated 
as a function of added salt concentration at fixed pH by computer simula- 
tions of the "primitive model" . The salt and counter-ions as well as a discrete 
charge pattern on the protein surface are explicitly incorporated. For param- 
eters roughly corresponding to lysozyme, we find that B2 first decreases with 
added salt concentration up to a threshold concentration, then increases to a 
maximum, and then decreases again upon further raising the ionic strength. 
Our studies demonstrate that the existence of a discrete charge pattern on the 
protein surface profoundly influences the effective interactions and that non- 
linear Poisson Boltzmann and Derjaguin-Landau-Verwey-Overbeek (DLVO) 
theory fail for large ionic strength. The observed non-monotonicity of B2 
is compared to experiments. Implications for protein crystallization are dis- 
cussed. 

PACS: 82.70.Dd, 61.20.Qg, 87.15.Aa 
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I. INTRODUCTION 



Interactions between proteins in aqueous solutions determine their collective behavior, in 
particular their aggregation, their complexation with other macromolecules, and ultimately 
their phase behavior, including phase separation, precipitation and crystallization. Any the- 
oretical analysis of the properties of protein solutions must rely on a reliable understanding 
of their interactions. A good example is provided by the control of protein crystallization, 
which is an essential prerequisite for the determination of protein structure by X-ray diffrac- 
tion Ijl],^ . While at present protein crystallization is still mostly achieved experimentally by 
"trial and error" , and on the basis of a number of empirical rules [§] , there is clearly a need 
for a more fundamental understanding of the mechanisms controlling protein crystallization, 
and this obviously requires a good knowledge of the forces between protein molecules in so- 
lution, and of their dependence on solution conditions, including pH and salt concentration 

sa-i- 

Protein interactions have various origins, and one may conveniently distinguish between 
direct and induced (or effective) contributions. Direct interactions include short-range repul- 
sive forces, which control steric excluded volume effects, reflecting the shape of the protein, 
van der Waals dispersion forces, and electrostatic forces associated with pH-dependent elec- 
tric charges and higher electrostatic multipoles carried by the protein residues [|^. Other, 
effective, interactions depend on the degree of coarse-graining in the statistical description 
and result from the tracing out of microscopic degrees of freedom associated with the solvent 
and added electrolyte, i.e. the water molecules and microions. Tracing out the solvent re- 
sults in hydrophobic attraction and hydration forces, while integrating over microion degrees 
of freedom leads to screened electrostatic interactions between residues, the range of which 
is controlled by the Debye screening length, and hence by electrolyte concentration. 

However, while coarse-graining through elimination of microscopic degrees of freedom, 
leading to state-dependent effective interactions is a priori a reasonable procedure to describe 
highly asymmetric colloidal systems, where particles have diameters of typically hundreds 
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of nm and carry thousands of elementary charges, this is obviously less justified for the 
much smaller and less charged proteins. In particular the assumption of uniformly charged 
colloid surfaces, leading to spherically symmetric, screened interactions between the elec- 
tric double layers around colloid particles, as epitomized by the classic DLVO (Derjaguin- 
Landau-Verwey-Overbeek) potential ceases to be a reasonable approximation at the 
level of nanometric proteins carrying typically of order 10 elementary charges. The reason 
is that length scales which are widely separated in colloidal assemblies, become comparable 
in protein solutions, while the discreteness of charge distributions on proteins can no longer 
be ignored, since the distance between two charged residues on the protein surface is no 
longer negligible compared to the protein diameter. Thus, electrostatic, as well as other 
(e.g. hydrophobic) interactions are much more specific in proteins, and must be associated 
with several interaction sites, rather than merely with the centers of mass as is the case for 
(spherical) colloidal particles. 

Another very important distinction between colloids and protein solutions is that the 
forces between the former may be measured directly, using optical means P|-pT]], while 
interactions between proteins can only be inferred indirectly, from measurements by static 
light scattering of the osmotic equation of state which, at sufficiently low concentrations, 
yields the second osmotic virial coefficient B2 , the main focus of the present paper. 

The variation of B2 with solution conditions yields valuable information on the underlying 
effective pair interactions between proteins. Moreover it was shown empirically by George 
and Wilson H that there is a strong correlation between the measured values of B2 and 



the range of solution conditions which favor protein crystallization [|T^,^,0. Only if the 
measured value of B2 falls within a well defined "slot" can crystallization be achieved. If 
B2 is too large, repulsive interactions predominate, leading to slow crystallization rates. On 
the other hand if B2 is highly negative, strong attractions lead to amorphous aggregation. 

The correlation between B2 and crystallization may be rationalized by noting that protein 
crystals generally coexist with a fairly dilute protein solution, the thermodynamic properties 
(and in particular the free energy) of which are essentially determined by B2. Coexistence 



between a dense solid phase and a dilute fluid phase is generally a signature of a very 



short-ranged attraction between particles as compared to their diameter [|T6|-|T^. 

For such short-ranged attractive interactions, the phase-separation into dilute and con- 
centrated proteins solutions expected on the basis of a mean-fleld, van der Waals theory, is 
in fact pre-empted by the freezing transition, i.e. the critical (or "cloud") point lies below 
the freezing line. The critical fluctuations associated with this metastable cloud point may 
lead to a signiflcant enhancement of the crystal nucleation rate [^, while the position of 
the cloud point in the concentration-temperature plane is strongly correlated with the virial 
coefficient B2 |[T6|| . 

The present paper focuses on the variation of B2 with ionic strength of added salt. 
This is a particularly important issue since "salting out " of protein solutions is one of 
the standard methods used to induce crystallization. An increase in salt concentration 
reduces the screening length and hence the electrostatic repulsion, allowing short range 
attractive forces (e.g. of hydrophobic or van der Waals origin) to come into play which 
will ultimately trigger nucleation. Recent experiments point to a non-monotonic variation 



of B2 with increasing ionic strength pT| , P^ , or to a pronounced shoulder in the B2 versus 



ionic strength curve [El] in lysozyme solutions. Closely related findings are the observation 



of a non-monotonic cloud point [p4H 26[|, and of a minimum in the solubility of lysozyme 



with increasing salt concentration ||2^; the solubility is obviously related to the osmotic 



virial coefficient |28|. Similarly, the attractive interaction parameter A, which controls the 



variation of the measured protein diffusion coefficient D with volume fraction, was found 



to exhibit a sharp minimum upon an increase of ionic strength of lysozyme solutions p9 
again, this interaction parameter strongly correlates with B2 

Traditional models for the protein-protein interaction cannot easily reproduce such non- 
monotonic behavior of B2 or related quantities. The "colloidal" approach based on spher- 
ical particles interacting via the screened Coulomb DLVO potential p] can only predict a 
monotonic decrease of B2 with ionic strength The same is true of models |]^,|T2|JT^ 



accounting for short-range attractions via Baxter's "adhesive sphere" representation |3^. In 



these models, which assume central pair-wise interactions, B2 reduces to a simple integral of 
the Mayer function associated with the spherically symmetric potential |3^J3^ . More recent 
calculations account for the asymmetric shape of proteins p^ , |36[] , or include several "sticky" 
sites at the surface of the protein p7| , |38| . 

In these traditional calculations, electrostatic interactions between proteins and mi- 
croions are routinely treated within mean-field Poisson-Boltzmann (PB) theory, generally in 
its linearized version (as is the case for the classic DLVO potential). However, as explained 
earlier, all relevant length scales (i.e. protein diameter, mean distance between charged sites 
on the protein surface, and between co and counterions, as well as the Debye screening 
length) are comparable in protein solutions, so that the discrete nature of both the inter- 
action sites, and of the co and counterions, can no longer be ignored. Moreover, Coulomb 
correlations are expected to be enhanced on protein length scales and may lead to strong 
deviations from the predictions of PB theory, which have recently been shown to induce 
short-range attractions, even between much larger colloidal particles p^p9| -|4^. 

The present paper takes into account the discrete nature of the microions within a 
"primitive model" description of the electrolyte, and presents results of Molecular Dynamics 
calculations of the equilibrium distribution of co and counterions around two proteins and 
of the resulting osmotic virial coefficient i?2- Two models of the charge distribution on the 
surface of the spherical proteins will be considered. In the colloid-like model the charge is 
assumed to be uniformly distributed over the surface, while in the discrete charge model, 
the charges are attached to a small number of interaction sites. The latter model will be 
shown to lead to a distinctly non-monotonic variation of B2 with ionic strength, as observed 
experimentally. 

The paper is organized as follows: The model and key physical quantities are introduced 
in section II. Simulation details are described in section III. Results of the simulations are 
presented and discussed in section IV, while conclusions are summarized in section V. A 
preliminary account of parts of the results was briefly reported elsewhere P5|. 
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II. MODELS, EFFECTIVE FORCES AND SECOND VIRIAL COEFFICIENT 



The globular proteins under consideration are modeled as hard spheres of diameter (T„, 



carrying a total (negative) charge —Ze. Within a "primitive model" representation [53 



the molecular granularity of the aqueous solvent is ignored, and replaced by a continuum 
of dielectric permittivity e, while the monovalent counterions and salt ions are assumed to 
have equal diameters Os and charges ±e. 

Two models are considered for the charge distribution on the surface of the protein. 
In the "smeared charge model" (SCM), the total charge —Ze is assumed to be uniformly 
distributed over the spherical surface, which is the standard model for charge-stabilized 
colloidal suspensions ||2^ , |39| -^, involving highly charged objects. According to Gauss' 
theorem, the SCM is equivalent to the assumption that the total charge Ze is placed at 
the center of the sphere. In the "discrete charge model" (DCM), point charges (— e) are 
distributed over a sphere of diameter cxrf = aOp (with a < 1, i.e. slightly inside the protein 
surface), in such a way as to minimize the electrostatic energy of the distribution. The 
resulting charge pattern, well known from the classic Thompson problem (see [^] for a 
recent review), is kept fixed throughout. Such Thompson patterns do not correspond to 
the true charge distribution on any specific protein (see where a simple toy model 

of lysozyme with different charge ditribution corresponding to solutions of different pH is 
constructed) but do provide a well defined discrete model for any value of Z. Note that the 
discrete distributions are characterized by non- vanishing multi-pole moments, depending on 
the symmetry of the distribution for any specific value of Z, while the SCM implies vanishing 
multipoles of all orders. 

At this stage the SCM and DCM models involve only excluded volume and bare Coulomb 
interactions (reduced by a factor 1/e to account for the solvent) between all particles, proteins 
as well as microions. 

The following physical quantities were systematically computed in the course of the MD 
simulations, to be described in the following section. 
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a) the density profiles of co and counterions around a single globular protein 

p±{r)=<J26{ff-f)> (1) 
j 

Here is the position of the j*'^ microion of species ± relative to the protein center, while 
the angular bracket denotes a canonical average over the microion configurations. For an 
isolated SCM protein these profiles are spherically symmetric, and depend only on the radial 
distance r = |r|. For isolated DCM proteins the profiles are no longer spherically symmetric, 
and may be expanded in spherical harmonics, as discussed in the Appendix. The anisotropy 
turns out to be weak, and only the spherically symmetric component (corresponding to 
averaging p±{r) over protein orientations) will be shown in the following. 

b) The second quantity, which will be the key input in the calculation of B2, is the 
microion averaged total force Fi = —F2 acting on the center of two proteins, placed at 
a relative position r = ri — the force Fi is a function of r. Its statistical definition 



was discussed earlier in the context of charged colloids |l39| , ^ , |49[| , and it involves three 
contributions: 

Fi = Ff) + Ff + (2) 
is the direct Coulomb repulsion between the charge distributions on the two proteins; 

~*(2) ~*(3) 

Fi is the microion induced electrostatic force, while Fi is the depletion force which may 
be traced back to the in-balance of the osmotic pressure of the microions acting on the 
opposite sides of protein 1 due to the presence of protein 2. F^ ' is directly expressible as 
the integral of the microion contact density over the surface of the protein . 

In the case of the SCM, the microion averaged force depends only on the distance r = \ fi2\ 
between the two proteins. For the DCM, on the other hand, Fi is a function of the relative 
orientations of the two proteins, as characterized by the sets of Euler angles Qi and ^2, i-e. 
F, = F,{f,Q,,Q2)- 

c) Once the force Fi has been determined as a function of r, Qi and ^2, one may then 
calculate an orientationally averaged, but distance resolved, effective protein-protein pair 
potential according to 



V{r) 



(3) 



where the angular brackets < ... >q^q^ refer to a canonical statistical average over mutual 
orientations of the two proteins weighted by the Boltzmann factor of the effective potential 
Veff{r,(li,fl2) such that dVeff{r,fli,0,2)/dr = —-^1(^5^15^2)- Explicitly, for any quantity 
A{f,Qi,Q2), 

_ J dnAAif, n^, ^2) exp{(-K//(r, n^, ^2)/kBT)} 



^ ^ ^^^'^^ J dnidn2 exp{-Veff{f, ni,n2)/kBT} 

d) The second virial coefficient B2 finally follows from the expression 

B2 = ^Jdf[l-b{r)] 



(4) 



(5) 



where 



6(r) 



dQidQ2 6xp(— Ve//(r, Qi, 0.2) /kBT). 



(6) 



The angular integrations are trivial in the case of the SCM, where V^ff depends on r. In 
the case of the DCM, one may use the identity 

d 



h{r) = exp 



(7) 



to show that B2 may be cast in a form similar to that appropriate for the SCM, namely 

B2 = \j df[l - exp{-y(r)/A;BT}] (8) 

where V{r) is the potential of the orientationally averaged projected force, as defined in 
Eq. @. As pointed out earlier, B2 is directly accessible experimentally by extrapolating 
light scattering data to small wavevectors ||35| or by taking derivatives of osmotic pressure 



data with respect to concentration [0,|Tj]. Results will be presented in the form of the 
reduced second virial coefficient B2 = B2/b'^^'' where b'^^'' = 27r(Tp/3, i.e. 

3 roc 



Q POO 

B; = 1 + — r^dr[l-exp{-Vir)/kBT}]. 



(9) 
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III. SIMULATION DETAILS 



We study a pair {Np = 2) of spherical proteins with center-to-center separation r, con- 
fined in a cubic box of length L = 4ap, which also contained monovalent co and counterions 
in numbers determined by their bulk concentrations and overall charge neutrality. There 
are ZNp counterions dissociated from the protein surface, and added A^^ salt ion pairs such 
that the screening of proteins is implemented via = Ng colons and = Ng + ZNp 
counterions in simulation box. A snapshot of a typical equilibrium microion configuration 
around two proteins is shown in Figure |I| for the protein charge number Z = 15. The two 
proteins were placed symmetrically with respect to the center along the body diagonal of a 
cubic simulation cell; periodic boundary conditions in three dimensions were adopted. L was 
chosen such that the box length is much larger than the range of the total (effective) protein- 
protein interaction, so that the results are independent of L for non-zero salt concentration. 
The long-range electrostatic interactions between two charged particles in the simulation 
box with periodic boundary conditions were modified using the Lekner summation method 



of images For our model to be a rough representation of lysozyme, we chose ap = Anm, 
and three different protein charges Z = 6, 10 and 15, corresponding to three different values 
of the solution pH. The microion diameter is fixed to be = Cp/lb = 0.267nm. 

For both the SCM and the DCM, the contact coupling parameter between a protein and 
a microion, namely F = 2e'^ /[ekBT{(Tp — ad + o-c)] for the DCM, and F = 2Ze'^ /[ekBT{ap + ac)] 
for the SCM, are comparable, and of the order of F ^ 3 at room temperature. We fixed 
the dielectric constant of water to be e = 81 and the system temperature to be T = 298i^r. 



Varying salt concentration for fixed protein charge Z corresponds to a fixed solution pH ^3 



Details of the runs corresponding to different salt concentrations are summarized in Table 
F Note that the Debye screening length td, defined by 



^87r(iV, + Z)g,e2' 
is less than lOA for salt concentration beyond O.IM. Here 
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V = V-'^{al + all (11) 

is the accessible volume for salt ions such that the salt concentration Cs is Ng/V' . Thus, 
the point charges on the protein surface are effectively screened from each other For 
each of the runs indicated in Table I, the distance-resolved effective forces and interaction 
potentials are calculated according to Eqs. (|2|,|D. The statistical averages over microion 
configurations leading to and Fi were evaluated from time averages in the Molecular 
Dynamics (MD) simulations. 



IV. MICROION DISTRIBUTIONS AROUND A SINGLE PROTEIN 

First, as a reference, consider a single protein {Np = 1) placed at the center of the 
simulation box. We calculated spherically averaged, radial microion density profiles p(r) = 
p+(r) + p-(r) in the immediate vicinity of the protein surface. For a single neutral sphere 
in a salted solution, results for p(r) are drawn in Figure |^. There is a marked depletion 
in the microion density, signaled by a minimum of p(r), well below the asymptotic bulk 
value. The depletion is enhanced upon increasing the salt concentration. At sufficiently 
high salt concentrations, this minimum is followed by a weak, but detectable, ion layer (see 
corresponding lines for runs 7 and 9 in Figure The formation of a depletion zone is not a 
consequence of the direct (hard core) interaction between salt ions and the protein surface, 
since the position of the observed layer is significantly further away from the protein surface 
than one ion diameter. A rough estimate for the distance between layer and neutral sphere 
gives a value of 2.5as, or equivalently O.lTcTp. For runs 7 and 9, where the ion layer emerges, 
this distance is of the order of an average ion separation in the system and twice the 
Debye screening length as well (see Table I). Obviously, it is the small ion correlations 
which lead to the peak formation in the salt density profiles. As an intuitive argument, 
we posit that the lack of mutual polarization in dense salt solution near neutral surfaces 
causes the ion depletion. Qualitatively similar depleted density profiles were observed in 
Lennard- Jones system confined between neutral planes ||5^ and in Yukawa mixtures 
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Furthermore, an effective force that pushes a single ion toward regions of higher sahnity is 



predicted within Debye-Hiickel theory for interfacial geometries pU 



Next we consider a protein sphere with charge number Z = 10. The density profiles 
of small ions are shown in Figure |^ for both the SCM and DCM. Figure ^ represents the 
corresponding total salt densities, as a sum of co- and counterion densities from Figure |^. 
At the lower salt concentration (up to run 5) the SCM and DCM models both yield an ac- 
cumulation of the micro ion density near ion-protein contact, in semi-quantitative agreement 
with the prediction of standard Poisson-Boltzmann (PB) theory. For rising ionic strength 
the total microion density gets depleted near the protein surfaces, as in the previously con- 
sidered case of a neutral sphere. Remarkably, this depletion occurs both with the SCM and 
DCM and contradicts the PB prediction. The intuitive picture is that a microscopic layer of 
counterions is formed around the proteins. An additional salt pair now profits more from the 
bulk polarization than from the protein surface polarization and is thus excluded from this 
layer. By normalizing the profiles to the total bulk density, this effect becomes visible as a 
depletion zone in Figure^, where a noticeable difference between the SCM and DCM profiles 
also emerges. Whereas the DCM predicts a contact value pc{r = {ap + ac)/2) larger than the 
bulk value, SCM predicts a much stronger microion depletion near contact. Together with 
this, the contact value of the DCM model is always larger than that of the SCM model for 
the same salt concentrations. This finding illustrates the sensitivity of correlation effects to 
the assumed charge pattern at the surface of a protein. This correlation effect is, of course, 
absent in the (non-linear) PB and DLVO theories, which always predict a monotonically 
decreasing density profile p(r) (see Figure]^). Within DLVO theory the density of plus and 
minus salt ions near protein surface in the SCM model are defined as 



Qs 47r r + p_ 



,^ ZpLvo kl e , 2p^p_ 
p_(r) = \ . (12) 

q, 471 r p^ + p_ 



where 
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N, + 2Z _ AT, _ exp(fcz3crp/2) 

P^-^^' ^^^^^ - ^TTW2"' 

kl = kl + kl, kl = ATcgyp^/eknT , kl = ATrgypJeksT . (13) 

Nonlinear PB equations were solved via iterations of the potential of a homogeneously 
charged sphere placed at the center of a spherical cell. The cell radius was determined 
by the given protein concentration. 

A direct comparison between SCM, DCM models and non-linear PB theory is shown 
in Fig. 1^ for two of the higher salt concentrations from Fig. ^. For the intermediate salt 
concentration Cg = 0.206 Mol/1 (run 4) both simulation and theory predict a monotonic 
decrease of salt density away from the protein surface, whereas, in the case of a dense salt, 
Cs = 0.824 Mol/1 (run 7), simulation results strongly deviate from the PB prediction. 

A multipole expansion of the total salt number density in the DCM, discussed in Ap- 
pendix A, demonstrates that the higher order expansion coefficients are strongly damped 
and much weaker than the zero-order homogeneous term shown in Figs. ^, § and ^. 

A. Effective force and B2 for a protein pair 

Next we calculate the angularly averaged effective interaction force F{r) = — ^^^^ and 
potential V{r) between two proteins embedded in a sea of small salt ions. Simulation results 
for the simpler case of the SCM, are plotted in Figure for Z = 10 and compared to the 
DLVO predictions. There is a systematic deviation between the theoretical and simulation 
results. While the DLVO theory ||^ potential 

^(i>Lyo)(^) = ^k^exp(-r/rz,), (14) 
er 

always results in repulsive forces, simulations indicate the possibility of an attraction between 
proteins for large salt concentrations. The force F{r) at the higher salt concentrations Cs, 
shows a maximum at a distance r nearly equal to the ion diameter. Note that, for the 
highest salt concentration considered, Cs = 2.061 Mol/1 (run 11), where the electrostatic 
interactions are almost completely screened out, the effective force F{r) is dominated by 
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entropic effects; it is reminiscent of entropic depletion force of hard sphere system. The 
corresponding potential is negative at short distances as shown in the inset of Figure 
and is related to the depletion in the microion total density profiles p(r) around an isolated 
protein, shown in Figure ^. We note that such an entropic attraction is not contained in 
DLVO theory. Its origin is also different from the salting out effect studied in P^J57| -pU[ or 



the macroion overcharging effect studied in |6T]|. In Figure | the salt dependence of the total 
interaction force F(r) (Eq.(H)) is broken down into its components F*-^-* and F^^^ for two 
values of r. This helps to show that at large salt concentrations, it is indeed the entropic 
component that causes the force to be attractive for run 11 in Figs ^ and ^ja. Finally, 
we mention that the range of attraction observed here will depend on the electrolyte (salt 
ion) size. This feature of our model may hint at a cause for the salt specificity observed 



in salting out experiments on protein crystallization [^]. In addition, as shown in ^3 
an entropic attraction from the electrolyte could lead even to phase separation in colloid- 
electrolyte mixtures. Including this term leads to an effective Hamaker constant - describing 
the dispersion interactions - that is lower, and in better accord with experimental findings 



The same calculations were carried out for the other two protein charges for the SCM 
model, Z = 6 and Z = 15, with qualitatively similar results to those obtained for Z = 10. 

It is clear that the effective forces and potentials between two proteins will no longer 
be spherically symmetric within the DCM model. For example, three distinguishable mu- 
tual orientations of the two proteins are schematically outlined in Figure ^a, corresponding 
to particular configurations of the Euler angles ^1,^2 of the two proteins. Nevertheless, 
our simulation results, presented in Figure Ipj, for the three orientations, show that the 
actual anisotropy of force is weak. However, this is no longer true for the Yukawa segment 
model P9| , |65| , |66[] , as shown in the inset of Figure ^o. Within this model, the total effective 
interaction potential between a pair of protein spheres is given by 

^^^'Hr) = 4 E ^^^"^°Hkl-r-;|), (15) 
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where Vk and r„ represent the positions of the point unit charges of different proteins. 
We emphasize that the aelotopic (or nonisotropic) interactions incorporated in our DCM 
differ from those considered, for example, in Ref. |^, where B2 is calculated for a set of 



hydrophobic attractive patches on protein surface. Within our version of the DCM the 
third configuration in Fig. ^ (solid line), has the highest statistical weight of the three 
cases pictured (see Eqs.(||,§)). If, on the other hand, the point charges on the protein are 
replaced by attractive patches [Q, then the configuration with two points nearly touching 



(dot-dashed line in Figure P), is the statistically most favorable conformation. Similar 
arguments hold within a molecular model for site-specific short-range attractive protein- 



protein interactions, [37,68|. Results for distance-resolved forces within the DCM model 



are shown in Figure |TO|a, for Z = 10. When the salt concentration is less than Cs ^ 0.2 



Mol/1, the results are similar to those of the SCM model: i.e. for low ionic strength, the 
force is repulsive, while for high ionic strength there is an attraction near contact followed 
by a repulsive barrier. The distinguishing property of the DCM is the nonmonotonicity of 
the force with the increase of ionic strength. This, in turn, gives rise to the nonmonotonic 
behavior of the spherically averaged interaction potential V{r) shown in Figure pUjb. This 
feature of V{r) manifests itself in the following way in Figure |TU|b: the potential is first 
strongly reduced as Cs is increased, then its amplitude and range increase very significantly 
at intermediate concentrations {Cs — 1 Mol/1), before it nearly vanishes at the highest salt 
concentrations. Note that V{r) even becomes slightly attractive at contact (r = dp) for 
Cs — 2 Mol/1. Similar effects are also observed for Z = 6 and Z = 15 (see Figure |lT]), 
suggesting that the effect is generic for discrete charge distributions. 

Once the effective potential V{r) is known, it is straightforward to calculate the second 
osmotic virial coefficient using Eq.(|D. In doing so, however, one should keep in mind that 
it is the total interaction that enters B2. Real proteins also exhibit an additional short- 
range interaction, as seen, for example, in experimental studies of the osmotic pressure and 
structural data for lysozyme or in fits to its phase-behavior |T^. This attraction stems 



from hydration forces, van-der-Waals interactions, and other molecular interactions that 
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are, to a first approximation, independent of salt concentration. Hence, we liave taken tlie 
expected sliort-range attraction between proteins into account by adding to tlie effective 
Coulomb potential in Eq.(P), an additional "sticky" sphere potential of the Baxter form 



VsHsjr) 



In 



12t(5 



ap<r <ap + 5 , (16) 



r>Op + 5 



with potential parameters 5 = 0.02crp and r = 0.12, which yield reasonable osmotic data 
for lysozyme solutions [|l^,^,^ in the high salt concentration regime. This square well 
potential is isotropic by nature and ignores the directionality in hydrophobic attraction 
between proteins [p8| , p7| . Short range attractions lead to "energetic fluid" behavior ffO 



where the crystallization is driven primarily by the details of the interactions, instead of being 
dominated by the usual entropic hard-core exclusions. This suggests that the directionality 
may be very important to details of the protein crystallization behavior ||3^. However, for 
the physically simpler behavior of the virial coefficient, the directionality can be ignored 
as a first approximation. For simplicity, we assume the parameter r to be independent of 
electrolyte conditions, although a weak dependence based on experimental observations is 



reported in ||23|j68|| . The addition of Vshs{^) strongly magnifies the nonmonotonicity of B2 
stemming from the nonmonotonic behavior of V{r) near contact. 

Results for B2 as a function of salt concentration are shown in Figure |12| for three different 
protein charges. There is a considerable qualitative difference between the predictions of the 
SCM and the DCM models for the variation of B2 with monovalent salt concentration 
Cs for each protein charge Z. Whereas the SCM (dashed curves in Figure [1^) predicts a 
monotonic decay of with Cs, the DCM leads to a markedly non-monotonic variation, 
involving an initial decay toward a minimum (salting-out) followed by a subsequent increase 
to a maximum (salting- in) and a final decrease at high Cg values (salting-out). The location 
of the local minima shifts to higher/lower values of Cs for larger /smaller protein charges 
Z. Thus for larger protein charge one needs a higher salt concentration to achieve the 
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"salting out" conditions conducive to protein crystallization . Even though the effective 



Coulomb potential between proteins is small, with an amplitude only a few percent of the 
thermal energy ksT, its effect on B2 is dramatically enhanced by the presence of the strong 
short-range attractive Baxter potential. 

We have also compared the effective potentials shown in figure O with a recently pro- 



posed scaling collapse of protein osmotic virial coefficients [^]. This scaling effect, observed 
for a number of experimental conditions , can be explained with simple arguments based 



on Donnan equilibrium |7l|. To lowest order, the effects of salt concentration and protein 



charge on B2 are to subsumed in the following approximate scaling relation: 

= 52 - Z^ACs, (17) 



where -62°^ is the bare virial coefficient, independent of charge effects. As shown, for ex- 



ample in Figure 1 of reference |j7T|, this simple relation hold remarkably well above a salt 
concentration of Cg ~ 0.25M for a wide range of experimental measurements of B2 for 
lysozyme, which all tend to a plateau value of -B^Z-Bg^'^^ ^ (—2.7 ±0.2). One implication of 
this observed scaling is that the attractive interactions that govern are indeed roughly 
independent of salt concentrations above Cs ~ 0.25M. When we applied the same scaling 



procedure to our B2 curves, a similar plateau develops for both the DCM and the SCM 
models, albeit with B^^ slightly less negative than that found in the experiments. One 
could, of course, very easily match to experiments by adjusting the value of r, but to keep 



contact with our earlier work ||4^, we don't do so. Clearly the scaling does bring the DCM 
and SCM ,82 's close together for a given Z, but for different Z (related to solution pH), 
the scaling collapse is not as good as that seen in experiments, since we observe a larger Cs 
before it sets in. Nevertheless, considering the high density of co and counter-ions in the 
simulation, it is remarkable that a simple Donnan theory based on ideal gas terms performs 
so well. 

The origin of the non-monotonic variation of B2 with Cs can be traced back to the 
subtle correlation effects which cause an enhancement of the effective Coulomb repulsion 
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at intermediate salt concentrations in the DCM. These effects cannot be rationahzed in 
terms of simple mean-field screening arguments [^. The protein-microion correlations are 
of a different nature to those in the SCM, where they lead to a much more conventional, 
monotonic decay of B2 with Cs, similar to that expected from a simple screening picture. 

In order to gain further insight into the physical mechanism responsible for the unusual 
variation of the effective Coulomb potential and of B2 with salt concentration in the DCM, we 
consider the influence of a second near-by protein on the microion distribution near protein- 
ion contact. We have computed the difference between "inner" and "outer" shell microion 



contact densities for Z = 10, as schematically illustrated in the inset to Figure |1J. The local 
microion density is no longer spherically symmetric, due to the interference of the electric 
double-layers associated with the two proteins. The difference Ap = — Pout between 
the mean number of microions within a fraction of a spherical shell of radius R = 0.6(Tp 
subtended by opposite 60° cones, is plotted in Figure |14| versus salt concentration. Ap is 
always positive, indicating that microions ( mainly counterions) tend to cluster in the region 
between the proteins, rather than on the opposite sides. This follows because they can lower 
the total electrostatic energy by being shared between two proteins. However, there is a 
very significant difference in the variation of Ap with salt concentration Cg, between the 
SCM and the DCM models. Both exhibit similar behavior for lower salt concentrations 
Cg < 0.5 Mol/1; for example, both show a small maximum around 0.2 Mol/1. But for salt 
concentrations above 0.5 Mol/1, the SCM predicts a monotonic decrease of Ap, while the 
DCM leads to a sharp peak in Ap for Cg — 1 Mol/1. This highly non- monotonic behavior 



clearly correlates with the non-monotonicity observed in Figs. |T0|, |TT] and |T^. The basic 
mechanism can be summarized as follows: For the DCM, the excess number of microions 
between the two proteins leads to an excess entropic pressure or force, as demonstrated in 
Figure |1^, which is the origin of the increased repulsion between proteins around Cg = 1 
Mol/1. The enhanced microion density arises from subtle crowded charge correlation effects 
that cannot easily be understood at a mean-field level. 
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V. CONCLUSION 



In conclusion, we have calculated the effective interactions and the second osmotic virial 
coefficient B2 of protein solutions incorporating the electrostatics within the "primitive" 
model of electrolytes. In this way we include nonlinear screening, overscreening, and corre- 
lation effects missed within the standard Poisson-Boltzmann description. For discrete charge 
distributions, the interactions and related B2 vary in a non-monotonic fashion for increasing 
ion strength while for the smeared charge model, a standard workhorse of colloidal physics, 
this effect is absent. These correlation induced effects are missed within non-linear PB the- 
ory, and similar coarse-graining techniques taken from the theory of colloids. In addition 
to this, our simulations indicate the necessity of taking entropic forces into account when 
treating systems on the nanoscale. These forces are believed to be essential in the salting-out 
effect |^2| and could lead to an attraction even between neutral globular proteins p4| . 

Our MD calculations can easily be extended to the more complex (pH dependent) charge 
patterns of realistic proteins [^. In fact, in some cases it may be easier to do a full MD 
simulation than to solve the non-linear PB equations in a very complicated geometry. We 
expect similar mechanisms to those found for the DCM to be active there, leading, for 
example, to an enhanced protein-protein repulsion at intermediate salt concentration. Since 
the second osmotic virial coefficient determines much of the excess (non-ideal) part of the 
chemical potential of semi-dilute protein solutions, we expect the non-monotonicity of B2 
to have a significant influence on protein crystallization from such solutions in the course 
of a "salting-out" process. The non-monotonic behavior also suggests the possibility of an 
inverse, "salting- in" effect, whereby a reduction of salt concentration may bring B2 into 
the "crystallization slot" 0,0 • The sensitivity of B2 to ion-correlation effects may help 
explain the salt specificity of the Hofmeister series |6^. Finally, we stress that our non- 



monotonicity is qualitatively different from that observed for added non-adsorbing |73|j74 
and absorbing polymers. 
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TABLES 

TABLE L Parameters used for the different simulation runs. Ns is tlie number of salt ion pairs 
in simulation box, Cs is the salt concentration in Mol/1, the Debye screening length is defined 



by Eqn.(|lOD and = y 47r(2^^+2Z) j average distance between salt ions for a given salt 

concentration. 



Run 


Ns 


Cs (Mol/1) 




rs/(rp 


1 














2 


125 


0.05 


0.34 


0.39 


3 


250 


0.103 


0.24 


0.31 


4 


500 


0.206 


0.17 


0.25 


5 


1000 


0.412 


0.12 


0.2 


6 


1500 


0.62 


0.1 


0.17 


7 


2000 


0.824 


0.085 


0.16 


8 


2500 


1.03 


0.077 


0.15 


9 


3000 


1.24 


0.07 


0.14 


10 


4000 


1.65 


0.06 


0.124 


11 


5000 


2.061 


0.054 


0.118 
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FIGURES 




FIG. 1. Snapshot of a typical MD-generated microion configuration around two proteins, 
separated by r = l.TcTp. Tlie proteins carry 15 discrete charges — e and the monovalent salt density 
is Cs = 0.206 Mol/1. The globular protein molecules are shown as two large gray spheres. The 
embedded small dark spheres on their surface mimic the discrete protein charges in the DCM 
model. The small gray spheres are counterions, while the black spheres are colons. 
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FIG. 2. Normalized total salt density profiles p{r) near single neutral sphere. pQ = Ng/V' is 
bulk density. The added salt concentration is increased from top to bottom (see the arrow which 
refers to the p{r) near protein surface) according to runs 1-5, 7, 9. 
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0.6 0.8 1 1.2 1.4 

FIG. 3. Rescaled density profiles of small ions near single protein surface for the SCM (a) and 
DCM (b) models. The protein charge is Z = 10. The added salt concentration is increased (shown 
by an arrow) from top to bottom for solid lines (counterions) and from bottom to top for dashed 
lines (colons), according to runs 2-5, 7, 9, 11. 
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FIG. 4. Same runs as Figure H but now for the total salt density near single protein surface. 
The arrow (a direction of added salt increase) applies to all runs except run 11, which is shown as 
a solid line with symbols. 
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FIG. 5. Same runs as in Figure ^, but now for DLVO theory (a), and non-linear PB theory (b); 
both are for the SCM. 
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FIG. 6. Total density profiles p{r) of salt ions around a single protein with Z = 10, for run 4 
(bottom set of curves) and run 7 (upper set of curves), comparing DCM simulations (solid line), 
SCM simulations (dashed line), and nonlinear Poisson-Boltzmann theory (squares connected by 
lines) . 



31 



-1.2 I ' ' ' ' ' 1 

1 1.1 1.2 1.3 1.4 1.5 1.6 




FIG. 7. Total interaction force F{r) (a) and interaction potential V{r) (b) versus dimensionless 
separation distance r/ap within the SCM for a protein charge Z = 10. The force is divided by 
Fq = ksT/Xs, where Xb = e^/eksT is the Bjerrum length. The added salt concentration is 
increased from top to bottom, according to runs 1-5, 7, 9, 11. Dashed lines correspond to the 
DLVO model. The inset in (b) shows in more detail the differences between the SCM simulations 
and the DLVO model potential for run 11. 
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FIG. 8. The total interaction force -F(circles) and its electrostatic, -F^^^ (squares) and en- 
tropic F^^^ (triangles) components versus salt concentration. The separation distance is fixed at 
(a) r/ap = 1 and (b) r/ap = 1.1. The simulations are for the SCM with Z = 10, and show that at 
high salt concentrations, the entropic force dominates. 
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FIG. 9. (a) An illustration of three different mutual orientations of two proteins. Points inside 
spheres represent protein charges in the DCM. (b) Total interaction force F{r) versus dimensionless 
separation distance r/ap for mutual orientations shown in (a) for run 5 and Z = 10 in the DCM. 
The inset shows the same, but for a Yukawa segment model. 
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FIG. 10. Total interaction force F{r) (a) and interaction potential V{r) (b) versus dimensionless 
separation distance r/ap for the DCM at Z = 10. Solid line- run 7, dashed line- run 8, dot-dashed 
line- run 9, open circles- run 11. The inset shows low salt concentrations, from top to bottom: runs 



1, 4, 5. 
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FIG. 11. The same as in Figure 10 but now for protein charge (a) Z = 6 and (b)Z = 15. The 
run numbers are placed next to corresponding curves. The result for run 1 is 4 times reduced in 
y-value to fit the y-axis scale. The inset in (b) shows low salt concentrations, from top to bottom, 
runs 1, 4, 5, 7. 
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FIG. 12. Normalized second virial coefficient B2 = 82/62^ of a protein solution versus salt 
concentration Cg- Results are shown for protein charges Z = & (dashed lines), Z = W (solid lines) 
and Z = 15 (dot-dashed lines). The lines with (without) symbols correspond to the SCM (DCM) 
model. Whereas the SCM virial coefficients decrease monotonically with increasing salt concen- 
tration, as expected from simple screening arguments, the DCM shows a marked nonmonotonic 
increase of B2 at intermediate salt concentrations. 
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FIG. 13. The same runs as in Figure 12, but now for the bare virial coefficient determined 



as — — (773^ — -. The arrow is a guide for eye for the direction of increasing protein charge Z. The 

^2 



scaling collapse at high Cg has been related to a Donnan equilibrium effect |^T|. Note that, on the 
scale shown, the nonmonotonicity is hardly visible. 
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FIG. 14. Difference in the microion density near contact Ap versus salt concentration for protein 
charge Z = 10 at a protein-protein separation of r = 1.2ap. The solid and dashed lines correspond 
to the DCM and SCM models respectively. The inset shows the angular range over which Ap 
is averaged (see text). The non-monotonic density profile for the DCM lies at the origin of the 
non-monotonic behavior seen for the forces, potentials, and virial coefficients calculated for this 
model. 
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FIG. 15. The protein-protein interaction force and its components at a protein-protein separa- 
tion r = 1.2(Tp, in units of ksT/XB, versus salt concentration Cg for a charge of Z = 10. From left 
to right: (a)- total interaction force, (6)- electrostatic component of interaction force, (c)- entropic 
component of interaction force. Solid line- DCM and dashed line- SCM results. This figure demon- 
strates that the difference between the two models arises primarily from the contributions of the 
entropic force. 
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APPENDIX A: MULTIPOLE EXPANSION OF COUNTERION DENSITY 

AROUND PROTEIN IN DCM 



We expand the numerically calculated counterion density around a protein with charge 
Z in a Laplace series of spherical harmonics: 

p(r, e,ip) = J2 Cnm{r)Pn (cosie)) cos{mip) + 5'„^(r)P™(cos(e)) sm{mip) (Al) 

n,m 

The multipole spherical expansion coefficients (MSEC) are calculated during an MD 
simulation via 



Coo 


i 






Cio 
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) cos(6'j 
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Cu 
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) sin(6'i 
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k) cos{6i) cos{ipi)), 
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3 sin(6 


h) cos{9i) sin(v9j)), 
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)3sin2( 


9i) cos{2ipi)), 


S22 




n 


)3sin2( 


9i) sin(2v9i)), 



where i runs over counterions. 

The protein charge is chosen to be Z = 12, so that the rotation symmetry axis through 
two surface charges has a fivefold symmetry. The xz and xy plane projections of the 
protein charge pattern are shown in Figure [T^ . 



The variation of the nonzero MSEC versus distance from the protein center are shown 



in Figure |17. 
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FIG. 16. A schematic diagram for the xy plane (a) and xz plane (b) projections of the protein 
surface charge distribution in DCM model for protein charge Z = 12. 




FIG. 17. Five nonzero multipole spherical expansion coefficients (MSEC) of counterion density 
p+{r) for the DCM with a protein charge Z = 12 and salt concentration Cg = 0.05 Mol/1. Note 
that the higher order expansion coefficients are magnified (xlOO). 
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